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Abstract 

This paper is devoted to the construction and analysis of a Moser- 
Steffensen iterative scheme. The method has quadratic convergence 
without evaluating any derivative nor inverse operator. We present a 
complete study of the order of convergence for systems of equations, 
hypotheses ensuring the local convergence and finally we focus our at¬ 
tention to its numerical behavior. The conclusion is that the method 
improves the applicability of both Newton and Steffensen methods hav¬ 
ing the same order of convergence. 

Keywords: Steffensen’s method, Moser’s strategy, recurrence relations, 
local convergence, numerical analysis 

Classification A.M.S. : 65J15, 47H17. 


1 



1 Introduction 

One of the most studied problems in numerical analysis is the solution of 
nonlinear equations 

F(x) = 0. (1) 

Many scientific and engineering problems can be written in the form of 
a nonlinear systems of equations, where F is a nonlinear operator defined 
on a non-empty open convex subset 0 of M m with values in M m . A powerful 
tool to solve these nonlinear systems of equations is by means of iterative 
methods. 

It is well-known that Newton’s method, 

{ xn given in O, 

’ _ ( 2 ) 
x n+ i = x n - F'{x n ) 1 F(x n ), n ^ 0, 

is the one of the most used iterative method to approximate the solution 
x* of (1). The quadratic convergence and the low operational cost of (2) 
ensure that Newton’s method has a good computational efficiency. But the 
existence of the operator F'{x n )^ 1 is needed at each step or equivalently 
to solve F\x n )(x n + 1 — x n ) = —F(x n ). It is also known another feature 
of Newton’s iteration: calculating inverse operators is not needed if the 
inverse of an operator is approximated; namely, to approximate P , where 
P E C(X,Y), the set of bounded linear operators from X into Y , we can 
phrase G(Q) = Q~ l — P = 0, so that Newton’s method is reduced to Q n +1 = 
2 Qn — QnPQrun ^ 0, provided that Q o is given. So, in na, Moser proposed, 
to solve nonlinear operator equations, the following iterative process: 

x 0 , B 0 given, 

< %n -\-1 — VCn B n F(^X n ^ : Tl P 0 

Bn+ 1 = 2 B n - B n F'(x n )B n , n > 0. 

This new iterative method, which can be considered as a Newton-type 
method, does not need to calculate inverse operators. It can be shown that 
the rate of convergence is -4^, provided the root is simple (see [lTJ). This 
process uses the same amount of information per step as Newton’s method, 
but it converges no faster than the secant method, so that this unsatisfactory 
from a numerical point of view. For that, Hald proposes in [16] the use of 
the following iteration: 
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( 3 ) 


XQ, Bo given, 

< *£n+l — X n B n F (^X n ^, Tl ^ 0 

B n+ 1 = 2R n - 5 n F'(x n+1 )B n , n > 0. 

Observe that the first equation is similar to the Newton’s method in 
which case B n is F\x n )~ l . The second equation is Newton’s method ap¬ 
plied to G{P) = -P -1 — F'(x n - |_i) = 0. We can stress two features of (4): it 
has the same rate of convergence of Newton’s method and it does not need 
to solve nonlinear equations at each step (it is “inversion free”). Moreover, 
from the convergence of the sequence x n , the convergence of the sequence 
B n can be given, so that it also produces successive approximations B n to 
the bounded right inverse of F' (x*) at the solution x*, which is very helpful 
when sensitivity of solutions to small perturbations are investigated. How¬ 
ever, as Newton’s method, this method has the same serious shortcoming: 
the derivative F'(x) has to be evaluated at each iteration. This makes it 
unapplicable to equations with nondifferentiable operators and in situations 
when evaluation of the derivative is too costly. The goal of this paper is to 
modify the previous two-step method in order to avoid the evaluation of any 
Frechet derivative. 

It is common to approximate derivatives by divided differences, so that 
iterative methods that use divided differences instead of derivatives are ob¬ 
tained. Remember that an operator [u, v\F\, u, v G H, is called a first order 
divided difference [20l I2T] , if it is a bounded linear operator such that 

[u, v; F] : fl C M m — > M m and [u, v; F](u — v) = F(u) — F(y). 

In this paper, we will start with Steffensen’s method HI 13 El El 0 El HE 
HH HU |19] to approximate a solution of the nonlinear system of equations 
F(x) = 0. It is a one-point iterative process given by the following algorithm: 

[ x Q given, 

< (4) 

( x n+ i =x n - [x n , x n + F(x n );F] L F(x n ), n > 0, 

The main interest of this iterative process lies in that the approximation of 
the derivative F'(x n ), that appears in each step of Newton’s method, by the 
divided difference [x n , x n +F(x n )- 1 7 ] -1 is good enough to keep the quadratic 
convergence |22j and, therefore, keeps the computational efficiency of New¬ 
ton’s method. However, in the Steffensen method remains an interesting 
problem to solve, its algorithm needs to solve a linear system of equations 
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at each step: 


| x_i ; x- 0 given, 

[^m *^n T P^n)i B] (*®n+l ®n) — B(x n fi fi. 0. 

In order to improve this fact, following the previous ideas applied in the 
case of Newton’s method, we consider the following Moser-Steffensen type 
method 


x 0 , B 0 given, 

< x n+ i =x n - B n F(x n ), n > 0 (6) 

Bn+i = 2 B n B n [x n -|_i, x n +i T F(x n +\)\ F^B n , 

where we have changed the resolution of the linear system in every step by 
making several matrix multiplications. 

Whereas the operational cost of both processes is similar, however the 
possibility of ill-conditioned linear systems appearing will be avoided simply 
taking appropriate initial matrix Bq , and therefore this previous algorithm 
will be able to be stable more easily. The matrix given by the divided 
difference [x n ,x n + F(x n );F], can be ill-conditioned and therefore the lin¬ 
ear system of equations previously indicated ([5]), in the classical Steffensen 
method Q, should cause instabilities. 

Remark 1.1 Along the paper, we consider the following first order divided 
difference inWL m : [u, w, F] = ([u, v; where 

[u,v;F]ij = —— (Ffim,... ,Uj,v j+ i,.. . ,v m ) - Ffim, ..., Uj_i, Vj, ... , v m )) , 

u j Vj 

(7) 

u = (ui,u 2 , ■ • .,Um) T and v = (yi,v 2 , ■ ■ .,v m ) T . 

The organization of the paper is the following. In Section 2 we study 
the local error equation, seeing his quadratic convergence. After, in Section 
3, we obtain a local convergence result for this iterative process. Moreover, 
in the Section 4, the numerical behavior is analyzed. From this study, we 
can conclude that the considered method ([6]) improves the applicability of 
Steffensen method 0. 

2 Local order of convergence 

We assume that F : Q C R m — y R m has at least 2-order derivatives with 
continuity on Q for any x € 11 lying in a neighborhood of a simple zero, 
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x* G fl, of the system F(x) = 0. We can apply Taylor’s formulae to F(x). 
By setting e k = x k — x*, the local order, and assuming that [F'(x*)] _1 
exists, we have 

F(x k ) = F(x* + ek ) = T(e k + 0{e \)) , (8) 

where T = F' (x *), e| = (e k ,e k ) G R m x R m . Moreover, from [T3] we 
obtain 

[^fc , 2 /fc ; T] = T (/ + A 2 {e k + e k ) + 02(efc, £&)), (9) 

where e k = y k — x* and A 2 = ^ T _1 F"(x*) G 2 2 (ft, R m )- We say that an 
operator depending on and e k is an 0 2 (e k ,s k ) if it is an O(e| 0 e^. 1 ) with 
go + 91 = 2, qi > 0, i = 0,1. 

If we expand in formal power series of e k and e k , the inverse of the 
divided difference given in ([ 9 ]) is: 

[x k ,y k \ F] _1 = (I - A 2 (e k + e k ) + 0 2 {e k ,E k ))T~ 1 . (10) 

These developments of the divided difference operator © were previously 
used in the study of Grau et a 1. m- In particular, for the operator Q k = 
[x k , x k + F{x k ) ; F] we take 

[x k , x k + F(x k ) ; F] = T (/ + A 2 (e k + (2 1 + T) e k ) + 0(e|)) , 

Note that we have the following relation between F(x k ) and Q k e k : 

F(x k ) = Q k e k + 0(el). (11) 

A first main local result is stated in the following theorem: 

Theorem 2.1 e k +i = E k e k + 0{e\), where E k = I — B k <d k . 

Proof. 

By subtracting x* from both sides of the first step of ([b]) and taking into 
account © we obtain 

£fc+i &k F k F^xjz') 

— B k Q k e k T 0(e k ) 

= E k e k + 0(e k ), (12) 

and the proof is completed. ■ 
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Other relations that are necessaries to obtain a second theorem are the 
following ones: 

Ek+i = I — B k+ 1 0fc_(_i 

= I — 2 B k 0fc + i + B k ©fc+i B k 0fc + i 
= (I-B k @ k+1 ) 2 . (13) 


I - B k 0fc + i = I — B k (0fc+i - Q k )B k 0 fc 
= E k - B k (Q k +i - 0fe) 

= E k — G k e k + e k ), (14) 

where G k = B k A 2 (21 + T) e £ 2(^5 R). 

A second main local result is stated in the following theorem: 

Theorem 2.2 E k+1 = (E k + G k e k ) 2 + 0 2 (E k , e k ). 

Proof. 

By substituting ( [14] ) into © we have 

E k +1 = (E k ~ G k e k ) + 0(E k e k , e|)) 

= E k + (G k e k ) 2 — E k G k e k — G k e k E k + o(E k ,E k e k ). (15) 

The proof is complete. ■ 


In a more precise way we can write Theorems 1 and 2 (see (12) and (15) 
respectively) in norm terms. Namely, 

llomll = 0(\\E k \\\\e k \\) (16) 

\\E k +i\\ = 0 2 (||^|U|e fc ||) = 0(||^|| 2 ,|| efc |||| 2 ,||T; fe |||| efe ||). (17) 

We have three possibilities: 


If ||£lfc+i|| = 0(\\E k \\ 2 ), then \\E k \\ = O (||g fc -i|| 2 ), and from © 
we take ||efc + i|| = O (||11 2 1|||), and applying (|16j) for k — 1, we 
have ||e fe || = O (||E fc _i||||e fe _i||). Hence, ||e fc+ i|| = O (plfe-i|| 3 ||efc-i||). 
Taking into account that from (16) with k — 1 we have ||efc|| 3 = 
O (||-E fc -i|| 3 ||e fc -i|| 3 ), finally, we obtain 

||e*+i|| = O (||efc|| 3 ||e fc _i||" 2 ) . (18) 
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The equation associated with (18) is p\(t) = t 2 — 3t + 2 = 0. The real 
positive root that coincides with the local order of convergence in this 
case is p = 2 (see [T8l 123] ). 


• If H-Efc+iH = O (||e fc |||! 2 ) and E k = O (||e fc — 1 1| || 2 ), then from (16) we 
have 

||efc+i|| = O (||efc|| ||efc—11| 2 ) • (19) 

In this case, the local order of convergence p is the unique real positive 
root of the indicial polynomial (see (MEHUS!) of the error difference 
equation (19) given by p(t) = t 2 — t — 2. That is, p = 2. 

• If H-Efc+iH = O (||£fc||||e fc ||) and we have (see ([16])) ||e fe+ i|| = O (\\E k \\\\e k \\), 
then 

\\ e k+i || = 0(\\E k+1 \\) = 0(||£; fc ||||^ fc _ 1 ||||e fc _ 1 ||) 

= O (||£'fc_ 1 || 2 ||eA ; - 1 || 2 ) 

= 0{\\e k f). 

The method, as in the two first cases, presents quadratic order of 
convergence again. 

We are now in a position to state the following theorem. 

Theorem 2.3 The iterative method defined by pi), from a local view of 
point, is a quadratic method. That is ||efc+i|| = O (|| e fcl| 2 )- M 


3 Local convergence analysis 

In this section, we prove the local convergence of the Moser-Steffensen type 
method given by ([6]). First, we establish a system of recurrence relations, 
from the real parameters that are introduced under some conditions for the 
pair (F,x o), where a sequence of positive real numbers is involved. After 
that, we can guarantee the semilocal convergence of the method in R m . 


3.1 Recurrence relations 

We suppose: 
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(Cl) Let x* a zero of nonlinear system F(x ) = 0 such that ||l ?/ (x*)|| < M 
and there exists r > 0 with B(x*,r) C Q. 

(C2) Let Bq G £(K m ,K m ) with Bq / 0 such that ||L>o|| = P and ||/ — 
B 0 F'(x*)\\ =6< 1. 

(C3) Consider xo G B(x*,r), such that there exists r > 0 such that xq F 
F(x o) G B(x*,r), with G B(x*,r) C J7. 

(C4) For each pair of distinct points x, y G 12, there exists a first order 
divided difference [x, y; F] of F and k ^ 0, such that 

llpy;^ - F'(x*)\\ < k(\\x- x*\\ F \\y -x*||); Vx,yGft, (20) 


Notice that, in these conditions, the Frechet derivative of F exists in 12 
and satisfies [x,x;F] = F'(x). On the other hand, to simplify the notation, 
we denote L n = [x n , x*; F] and 0 n = [x n , x n + F(x n )\ F], for all n G N. 

From the above, we denote «o = r, Po = P, So = S, do = r and define 
the scalar sequences: 


— (S n — i F k (Sn—iCXn—i^CXn—i, Ot n — (1 “I - M F k (y n )CY n 
dn — S n F k fd n [oi n j r \ F dn+i) ; Pn — (1 "F d n —\)f3 n —\ (21) 

Sn = Sl_ 1 + k MPn_l{ a n + dn) 


Next, for n 
pTj) and {x n }: 


1 , we prove he following recurrence relations for sequences 


\\xi - x*\\ 

< 

a l, 

F(x i)-x*|| 

< 

di, 

||/-5,0211 

< 

di, 

Pill 

< 

Pi: 

~BiF\x*)\\ 

< 

Si, 


provided that 


x 0 ,x 0 F F(x q) G 12. 
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From the initial hypotheses, it follows that x\ = xq — BqF(xo) is well 
defined and 

||*i - x*\\ = ||x 0 - B 0 F(x o) - x*|| 

= ||x 0 - X* - B 0 L 0 (x 0 - x*)|| 

< ||/-S 0 Lo||||*o-*l 

< (II I- Bo *V)II + PollULo - F'(x*)ll) ||x 0 - *1 

< (A) + A) ka 0 ) ||x 0 - x*|| 

< (A) + A) ka o) a 0 = ai. 


On the other hand, 

||.Ti + F(x i) - x*|| = ||xi + Li (xi - x*) — x*|| < II/ + Li||||xi - x* 

<(||/ + F'(x*)|| + imx*)-i 1 ||)||x 1 -xl 

< (1 + ll^ ?/ (x*)|| + k || xr — x* ||) ||xr — x* || 

< (1 + M + k || xi — x*||) ||xi — x* || 

< (1 + M + k ot\) ot\ = oi\. 


Assuming that a\ < ao and (1 + M + kr)r < f, then a\ < a o and 
therefore xi,xi + F(xi) £ 0. So, there exist ©i = [xi,xi + F(xi)\ F] and 
B\ = 2Bq — BqQiBq. Then, we establish 

||/ - Bo©!|| < 11 / - B 0 F'(x*) II + ||B O ||||01 - F l (x*)\\ 

< ||/ - B 0 F'(x*)\\ + H-Boll|| fc(||x! - x*|| + ||xr + F{x\) - x*||) 

< So + /do k(a i + ai) = d 0 . 


As a consequence, 

||Bi|| = H2S0-BoQiBoll 

< (1 +l|/-B 0 © 1 ||) ||B 0 || 

< (l + do)||B 0 || = An 
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Moreover, to finish the first step, notice that 


||/ - B 1 F'(x *)|| = ||I - ( 2B 0 - B 0 QiB 0 ) F'(x*)\\ 

< \\I- B 0 F'(x*)W 2 + ||5o|| II©! -F\x*)\\ \\B 0 F\x*)\\ 

< ||/ - B 0 F\x*)\\ 2 + pof ||F'(x*)ll k (||a:i - x*|| + H*! + F( X1 ) - **||) 

< + 0o M k (ai + ai) = p 


Next, for n 
(21) and {x n }: 


2, we prove the following recurrence relations for sequences 


11*2 -*1 

< 

a 2 , 

||*2 + F{x 2 ) - **|| 

< «2, 

||/-F 2 0 3 || 

< 

d 2 , 

\\b 2 \\ 

< 

P 2 , 

\\I-B 2 F\x*)\\ 

< 

s 2 . 


Now, by ((gJ) , X 2 = x\ — BiF(xi), furthermore 

\\X2 - X*\\ < ||I - Bi [*!,**; F]||||*1 - **|| 

= ||Z — B\ L \||||*i — x* || 

< (||/ - B 1 F'(x*)\\ + Pillpi - F'(x *)||) ||*i - x* 

< (<5i + j3\ k cti) ai = a 2 . 


Moreover, it is easy to check 

11*2 + F(x 2 ) — **II = 11*2 + B 2 (x 2 - X*) — ** II < ||F + F>2II 11*2 - X* 

<(I|/ + f'P)|| + ||f'(x*)-l 2 ||)||x 2 -x*|| 

< (1 + ||F'(**)|| + k \\x 2 - **||) ||* 2 - **|| 

< (1 + M + k ||* 2 - **||) ||* 2 - ** || 

< (1 + M + k a 2 ) a 2 = d 2 . 


Assuming that a 2 < «i, then d 2 < a\ and therefore x 2 ,x 2 + F(x 2 ) G fh 
So, 0 2 and B 2 are well defined. 
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Then, we have 

/-£10211 < ||i-£i£V)II + I|£ 1 |II|02-£V)II 

< ||I - £i£ / (x*)|| + IlSillll fc (||® 2 - x*|| + ||x 2 + F(x 2 ) - x *||) 

< (5i + /?! k(a 2 + a 2 ) = d±, 


and we get 

||£ 2 || = ||2 j b 1 -b 1 0 2j b 1 || 

<(l + ||/-£ 1 0 2 ||)||£ 1 || 
<(l + di)||Bi|| 

< (1 + di) Pi = p 2 . 


To hnish the second step, we consider 
I-B 2 F\x*)\\ = \\I-(2B 1 -B 1 e 2 B 1 ) TV)|| 

< II* - £i F '( x *)lf + \\Bi\\ ||02 - £'(**)ll H-Bi T'(x*)ll 

< ||/ - Bi T'(x*) || 2 + H-Br|| 2 ||i ?, (x*)|| k (\\x 2 - x*|| + \\x 2 + F(x 2 ) - x*||) 

< Si 2 + Pi 2 M k (a 2 + a 2 ) = 5 2 . 


At this time, we are able to obtain a general result that allows us to 
relate the sequences (21) and {x n }. 


Lemma 3.1 In the previous conditions, if the sequence {a n } is decreasing 
and (1 + M + kr)r < r, then 


(I) 

\\x n 

- x* < a n , 

(II) 

| X n 

+ F(x n ) - x* | < a n , 

(III) 

Wi¬ 

B n @n+l\\ P d n , 

(IV) 

ll Bn 

II ^ Pn> 

(V) 

II /- 

-B n F’(x*)\\<6 n , 


for all n e N. 
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Proof. 

First, since {a n } is a decreasing sequence and (1 + M + kr)r < r, then 
{an} is also decreasing. Therefore, we have that Xk,Xk + F{xk) G fi is 
verified for k = 0, 1 ,..., n, and therefore there exist ©*,, Bf. such that 
Xk+i = Xk — BkF(xk) is well defined. 

Once shown the relationships (I) — ( V ) in their first two steps previously, 
we will use mathematical induction. Suppose that the relations (/) — (V) 
are true for k = 1,..., n and we are going to show them for k = n + 1. 

Observe that 


F{x n ) = F{x n ) - F(x*) = L n (. x n - X*) 

and by Q 

Xn +1 X — X n B n F (r n ) X — (/ B n LnfiXn x ), 

furthermore 

||^n+i x || ^ ||I B n L n || ||x n x || 

< (\\I-B n F'(x*)\\ + \\B n \\\\L n -F'(x*)\\) \\x n -x*\\ 

< (||/ - B n F'(x*)\\ + ||5 n || k \\x n - x *||) ||x n - x*\\ < ( 5 n + (3 n ka n ) \\x n - x*|| 

is (fin T fin k On) — C^ra+l- 

In addition 

\\x n +i + F(x n+ 1 ) - x* || < \\I + L n+ 1 1| ||x n+ i - x* || 

< {\\I + F'(x*)\\ + \\F'(x*)-L n+1 \\) \\x n+1 -x*\\ 

< (l + ||F'(x*)|| + k \\x n+1 - x* ||) ||x n+ i - x* || 

< (1 + M + k ||x n+ i - a;*||) \\x n+1 - x*|| 

< (1 + M + k a n +i) a n +i 
= <5n+ 1 • 

Assuming that a n+ \ < a n , then d n+ i < a n and x n+ \, x n+ \ + F(x n+ 1 ) G 
fh Consequently 0 n +i and B n+ \ are well defined. 

So, we can consider 

||/ - -BnOn+ill < IIJ - B n F\x*) + ||B n ||||0 n+1 - F'(x*)\\ 

< \\I - B n F\x *)|| + ||B n |||| k(\\x n+l - x*|| + \\x n+1 + F(x n+ 1 ) - x* 
is fin T Pn k(a n j r 1 T d n +i) — d n , 
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that implies 

H-Bn+lH = ||2S n -S n 0 n+ iS n || < (1 + ||/ — -B n © n _|_i||) \\B n \\ < (1 + d n ) fi n = fin+1- 

Finally, to round off 

||/ - B n+1 F'(x*)\\ = ||I - (2 B n - B n @ n+1 B n ) F'(x*)\\ 

< ||/ - B n F'(x*)\\ 2 + \\B n \\ ||0 n+1 - F'(x*)\\ \\B n F\x*)\\ 

< ||/ - B n F'(x*)\\ 2 +\\B n \\ 2 ||F'(^)|| k (Hsn+i - x*\\ + ||x n+ i + F(x n+1 ) - x*\\) 

< ^rF + fin M k (a n +i + d n -|-i) = 

The proof is complete. ■ 

Once generalized the previous recurrence relations to every point of the 
sequence {x n }, we have to guarantee that {x n } is a Cauchy sequence having 
into account these recurrence relations. For this, we first analyse the scalar 
sequences given by(|21|) in the next section. 


3.2 Analysis of the scalar sequence 


Now, we analyse the scalar sequence defined in (21) in order to prove later 
the semilocal convergence of the sequence {x n } in M m . For this, it suffices 
to see that {x n } is a Cauchy sequence. First, we give a technical lemma. 


Lemma 3.2 Let {a n }, {a„}, {d n } and {<5 n } be the sequences given by 
pip . If it is verified that 

(1 + M + kr)r < r, hi < 5o and (1 + do) 2 (^o + kfiafi) < 1, (22) 

then 

(a) (h 0 + kfi 0 a 0 ) < 1 and (1 + doX^o + kfi 0 a 0 ) < 1, 

(b) the sequences {a n }, {n}, {d n } and {h n } are decreasing. 


Proof. 

Observe that as (1 + do) > 1, then (1 + do) 2 (S + kfia) < 1 implies that 
(1 + do)(5 + kfia) < 1. By the same reason (h + kfia) < 1 and (a) holds. 
We shall prove (b) by induction. 

From (a), for n = 1, we have that a\ = (ho + k fio «o) “o < <^0 and 
di = (1 + M + k an) ai < (1 + M + k ao) ao = ao- 
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For n = 2 , having into account that P\a.\ < (l + do)A) (do + k A) a o ) a o < 
A)Oo, we obtain «2 = (A + A Al a\) a\ < (< 5 o + k A) ao) ao = ai and 02 = 
(1 + M + k 02 ) a2 < (1 + M + k a\) a\ = a\. 

To analyze the sequences {d n } and {d n }, we must also have in mind that 

/0i«2 < (1 + d 0 ) 2 /3o0 2 < /3o(l + do) 2 (A + A A. ai) ai 

< 0q{\ + d 0 ) 2 (d 0 + kp 0 a 0 )ai < $,01 

and therefore it follows that /3 2 a2 < / 3 gdi, A1O2 < A)«i and A1A2 < A)Oi- 
Then, to finish the case n = 2 , taking into account that by hypothesis 
< 5 i < 5 q , we get 62 = S 2 + k Mfi\{a.2 + 0:2) < < 5 $ + k Mp^{a± + di) = A and 
di = 5 i + k A( a 2 + d2) < d + k/ 3 (ai + di) = do- 

From now, we suppose that ao > ai > ■ ■ ■ > a n , A) a o > Al q i > ■ • • > 
/3 n a n and fi^ai > P 2 ot 2 > • • • > /3 2 a n +1 hold, which implies that the se¬ 
quences {d fc }£ =1 , {/3fcd fe }^ =1 , {Pi a k+ i}fc =1 , {/3fcafc+i}fc =1 and {Afcd fc +i}£ =1 
are decreasing, as well as {d k } k= 1 and {dfc}^~g are. 

We need to check the inductive step. 

In first place, it is easy to prove: 

a n +l — {dn T k P n a n )a n <C {$n —1 T k Ai— l&n— l)Ojj— 1 — O n , 
a n +i — (1 T M -)- k (x n -\-i )on+i <1 (1 T AI" -)- Aa n )an — a^, 
dn = d n + A Ai(On+l T dn+l) <1 d n — \ + A Ai —l(On T dn) = d n —h 
Ai+i = d 2 + A: Af/3^(a n +i + d n +i) < d 2 _! + A’ M /3 2 _i(a n + d n ) = S n . 

On the other hand 

Ai+iOn+i ^ A(1 + d n )a n+ i < a n -w- (1 + d n ){d n T A/3 n a n )a n <C a n 

(1 + d n ){6 n + A AiOn) ^ 1) 

which is true since that 

(1 + d n ){d n + A’ Pn a n) < (1 + do) (5 + k /3 qlq) < 1. 

Note that we also have 

Oji+2 — (A 1+1 T A/3n+io n +i)o n +i <C (d n + k f3 n a n )a n — o n +i 

and therefore, 

Ai+ 1 0++2 <C (1 T d n ) P n a n +2 ^ AA 1 ~f~ d, / ) (dn-f-1 T A 3n -\-1 O n+ i) a n +i 
^ Aj(l T d n ) {d n T k (3 n ot n ) a n <C AiOn+i- 

Consequently {o n }, {d n }, {d n } and {d n } are decreasing. 

The proof is complete. ■ 
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3.3 A local convergence result 


First of all, we notice that, fixed xq G B(x*,r), if there is r, with the 
condition (1 + M + kr)r < f (in order to verify B(x*,r ) C fi) and, in 
addition, r verifies the two other conditions given in ( 22 ), then {x n } will be 
convergent. 

On the other hand, the conditions given in (22) can be written in the 
following way: 


(1 + M + kr)r < f, (23) 

<5i < <5 q d 2 + kMf3 2 (ot\ + di) < (5 0 < <5(1 — 5) — rPi(r), (24) 


(1 + do) 2 (^o + k/3ao ) < 1 (1 + 6 + kf3{a\ + di)) 2 (<5 + k/3r) < 1 

<^0< (l-(l + <5) 2 5)-rP 3 (r), (25) 

where Pi is a polynomial of degree one and P 3 is a polynomial of degree 
three, both decreasing and concave in (0, + 00 ). In this situation, it is clear 
that if 1 — (1 + 5) 2 5 > 0, then there is always r checking (24) and (25). 


Theorem 3.3 

Using the above notations, under the initial conditions (Cl) — (C4), we 
assume that there exists r > 0 verifying the conditions (25 1 ) and 

B(x*,r) C 12. Then, if we consider xq G B(x*,r), the sequence { x n } given 
by (ffy is well defined and converges to a solution x* of F(x) = 0. 


Proof. 

First, it is easy to prove from the hypothesis that x n ,x n + F(x n ) G 0 
for n > 1. Then, the sequence {x n } given by © is well defined. 

On the other hand, if we denote L = 5 + k/3r , we have 


Haq — x*|| < (<5 + k/3r)\\xo — x*|| = L\\xo — P|| 

\\x 2 — a;*|| < (di + kfiioti)\\xi — x*|| < (5 + kf3r)\\xi — x*\\ < L 2 \\xq — x* 
So, in general, we obtain 


||x n -x*|| < (5 n -i+k/3 n -ia n -i)\\x n -i-x*\\ < (5+kj3r)\\x n -i~x*\\ < L n ||x 0 -a:* 

Since L < 1 from lemma 3.2, it follows that the sequence {x n } given by © 
converges to a solution x* of F(x) = 0. ■ 


Notice that, if O = M m , the condition (23) it is not necessary. 
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3.4 An example 

Next, we illustrate the previous result with the following example given 
in [10] . We choose the max-norm. 

Let F : M 3 —> M 3 be defined as F(x,y,z) = ( x,y 2 + y,e z — 1). It is 
obvious that the unique solution of the system is x* = (0, 0, 0). 

From F, having into account 0 , we have 

/ 1 0 0 \ 

F\x,y,z) = 0 2y + l 0 

\ 0 0 e z ) 

( l ° 

[(x,y,z),(u,v,w)-F]= \ 0 y + v+l 

\ 0 0 

So, F'(x*) is the identity matrix 3x3. Then, ||.F / (a:*)|| = 1 and consider 
Bq = diag{/3, j3, /?} with (3 = 0.75 and <5 = 0.25. On the other hand, there 
exists r = 1, such that B(0,r) = {tc G M 3 : ||u>|| < 1} C M 3 , and it is easy 
to prove that 

pZ _ e w 

\\[(x,y,z),(u,v,w)-F] -F'(0,0,0)|| < max{|y + u|, |----1|} 

z — w 

< \\{x,y,z))\\ + ||(u,u,ui)||, 

in B(0,r). Therefore, M = 1, k = 1 and considering ao = 0.246627 we 
obtain: 

(1 + M + kr)r = 0.554078 < f, = 0.107275, di = 0.226058, 

= 5.55112 x 10“ 17 , c7 0 = 0.5, 1 - (1 + d 0 ) 2 (5 0 + k/3a 0 ) = 0.0213177 

Therefore the iterative process 0 is convergent from any starting point 
belonging to B(x*, 0.246627). 

4 Numerical experiments 

In this section, we include two experiments to test the proposed algorithm 
([6]). In the first one, we check numerically its order of convergence and its 
stability behavior. Moreover, we propose specific chooses of the initial ma¬ 
trix Bq in order to improve the stability of the classical Steffensen method 
Q. In the final example, we apply our iterative method to solve the non¬ 
linear systems of equations that appear to approximate a stiff differential 
problem with an implicit method. 
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4.1 Academic example 

We consider the academic system of equations F(x, y ) = (0, 0) given by: 


X ^ np‘ 

(2x--) + (y- y -) = 0, 

x + y = 0. 


(26) 


This system has as solution (x*,y*) = (0,0) and its Jacobian matrix is 
given by the following matrix 


2 — 2- 1-K 

e e 


On the other hand, having into account 0, the matrix [x,x;F\, with 
x = (xi,yi) and x = (x 2 , 2 / 2 ), is given by 

/ 9 — x l+ x 2 1 _ V2 J ry2 \ 

( l' "l‘ ) 

• Stability 

The parameter e (when e —>• 0) increases the condition number of the 
difference divided matrix for a given initial guess and the Steffensen 
method Q should have problems of convergence to the solution for 
small values of e. Moreover, notice that F (e, e) is not invertible. 

We compute the maximum of the condition numbers of the linear 
systems that appear in the application of the Steffensen method Q 
(||A|| • ||A _1 ||) and the maximum of the conditions numbers in all 
the matrix multiplications in the Moser-Steffensen type method (|6| 

1 \\ab\\ >■ 

In table 1, the vector (e, e) is inside of the ball containing the initial 
guess and the solution. The Steffensen method Q diverges. 

This numerical behavior is similar for smaller parameters of e, as we 
can see in the table 2 (e = lO^ 1 ). The maximum of the condition 
numbers for the Steffensen metho d (|4[ ) is 6.09 10 2 (too big) and for 
the Moser-Steffensen type method ([6]) smaller than 30. Moreover, the 
sequences of condition numbers for both methods are increasing and 
decreasing sequences. 

For this example, the Steffensen method 0 has a small region of 
convergence. However, the method ([6]) only reduces a little its velocity. 
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Table 1: System (26). Errors for (xo,yo) = (—1,1), 1,1)|| < 10 3 

and e = 1. 


Table 2: 

|| B 0 F '(-0 


n 

Steffensen 

Moser-Steffensen 

i 

1.55 10 1 


1.00 10° 

2 

1.84 10 1 


3.76 10 1 

3 

2.12 10 1 


1.08 lO” 1 

4 

2.41 10 1 


1.80 10 -2 

5 

2.69 10 1 


9.12 10" 4 

6 

2.98 10 1 


3.61 HP 0 

7 

3.26 10 1 


7.60 10 -11 

8 

3.54 10 1 


4.18 10" 2U 

n (26). Errors 

< ICE 3 and e = 10' 

"o 

o 

A 

fH 

£ 

1 

n 

Steffensen 

Moser-Steffensen 

1 

1.70 10° 


6.84 10" 1 

2 

1.98 10 u 


2.08 10" 1 

3 

2.27 10 u 


8.90 10“ 2 

4 

2.55 10 u 


3.95 10" 2 

5 

2.83 10 u 


6.80 10 -2 

6 

3.12 10 u 


6.95 10" 4 

7 

3.41 10° 


1.22 10 -5 

8 

3.69 10 u 


5.38 10 -9 

9 

3.97 10° 


1.33 10 -15 

10 

4.25 10 u 


1.00 10” 28 


(-0.25,0.25), 


This is the main advantage of the Moser-Steffensen type method Q. 
The condition number of the operations used in the classical Steffensen 
method Q should be large (in the previous cases of divergence the con¬ 
dition number goes to infinity) and there is not a general strategy to 
find preconditioners for a given linear system. However, by construc¬ 
tion, the condition number in the operations of the proposed method 
([6]) (matrix multiplications) seems controlled with our election of Bq. 
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con- 


In Moser-type algorithms, as ([hj, the sequence of matrices B n 
verges to the inverse of the Jacobian at the solution. For this reason, 
a good candidate for Bq is an approximation to the inverse of the 
Jacobian at the initial guess [12]. 

• Order of convergence 

In table 3, the method obtains the results expected by our theoretical 
analysis and the second order convergence is clear. The vector (e, e) 
is outside the convergence region associated to the initial guess. We 
compute the maximum of the condition numbers of the linear systems 
that appear in the application of the Steffensen method Q and the 
maximum of the conditions numbers in all the matrix multiplications 
in our method ([b]). In this case, the maximum for the Steffensen 
method @ and for our method ®> are smaller than 10, and both 
methods work well. 


Table 3: System (26). Errors for (x'o, yo) = (—1,1), 
and e = 3. 


\B 0 F'(-1,1)\\ < 10 


-3 


n 

Steffensen 

Moser-Steffensen 

1 

1.41 10° 

3.55 10" 1 

2 

7.61 10" 1 

5.09 10“ 2 

3 

2.40 10" 1 

1.86 nr :i 

4 

2.60 10~ 2 

3.74 10“ e 

5 

3.11 10” 4 

2.01 10~ n 

6 

4.56 10 -8 

7.10 10~ 22 

7 

9.79 IQ” 16 



The results in the table 4 are similar to the first case, both methods 
have second order of convergence. 

In table 4 case, we consider a smaller e but again the point (e, e) is 
outside the ball including the initial guess (closer also to the solution) 
and the solution. In this case, our method preserve the second order 
of convergence (condition number smaller than 10), however the Stef¬ 
fensen method diverges (with condition number 1.80 10 2 in the last 
iteration). 

• Election of the initial matrix Bq 
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Table 4: System (26). Errors for (. xo,yo ) = (—0.5,0.5), 0.5,0.5)11 < 

10 -3 and e = 1. 


n 

Steffensen 

Moser-Steffensen 

1 

1.49 10° 

2.12 10 -1 

2 

2.75 10° 

4.21 10 -2 

3 

3.35 10° 

3.09 10 -3 

4 

6.82 10° 

2.68 10 -5 

5 

9.80 10° 

2.77 10“ 9 

6 

1.27 10 1 

3.75 10- 17 

7 

1.56 10 1 

8.71 10" 33 

8 

1.84 10 1 


9 

2.13 10 1 


10 

2.41 10 1 



In general, as we indicated before, a good candidate for Bq is an ap¬ 
proximation to the inverse of the Jacobian matrix at the initial guess. 

On the other hand, it is not necessary to consider Bq a really accurate 
approximation to the Jacobian at the initial guess, as we can see in 
table 5. In particular, we can take some iterations of some of the 
algorithms that appear for instance in (2]. 


Table 5: System (26). 
e = 3. 

Errors for Moser-Steffensen, 

{xo,yo) = (-2,2) and 

n 

|T?oF'(-2,2)|| < 1 \\BqF'(—2,2)\\ < 10 1 | 

|£? 0 F'(-2,2)|| < 10- 3 

i 

1.84 10° 

1.10 10° 

9.44 10~ 4 

2 

8.75 10' 1 

2.83 10 1 

2.24 10 1 

3 

2.53 10 i 

3.56 10 -2 

2.43 10 -2 

4 

3.29 10 -2 

9.86 10 -4 

4.87 10 -4 

5 

8.61 10” 4 

1.12 10 -6 

2.82 10" 7 

6 

8.37 10" 7 

1.89 10” 12 

1.21 10~ 13 

7 

1.01 10” 12 

6.57 10 -24 

2.71 10' 2ti 

7 

1.83 10“ 24 




Finally, in the table 7, we force the method to take the bad iteration 
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(e, e) (here the Jacobian is not invertible). Only with the new approach 
we are able to find the solution (the Steffensen method diverges). 

Indeed, we consider (xo> 2/o) = (e, e). For this situations, a possibility is 
to take the initial matrix Bq = 5I 2 , where 5 is a small parameter (10 -2 
in our experiments) and I 2 is the identity matrix. When the method 
leaves the conflict zone, it recovers its good properties (second order 
convergence). 


Table 6: System (26). Errors for (zo,2/o) = (2,2), e 


n 

Moser-Steffensen 

10 

1.13 10 2 

11 

2.81 10“ 4 

12 

2.07 10~ 7 

13 

1.30 10 -13 

14 

5.88 lO" 26 


2 and Bq 


IQ " 2 I 2 . 


4.2 A stiff problem: Chapman atmosphere 

This model represents the Chapman mechanism for the generation of the 
ozone and the oxygen singlet. In this example, the concentration of the 
oxygen y 3 = [O2] will be held constant. It is a severe test for a stiff ODE 
package [15] governed by the following equations: 

y[(t) = 2 k 3 (t)y 3 + k A (t)y 2 (t) - (kiy 3 + k 2 y 2 (t))yi(t ), 

y' 2 (t) = kiyi(t)ys - (k 2 yi(t) + k A (t))y- 2 (t ), 

with y 3 = 3.7 x 10 16 , h = 1.63 x 10“ 16 , k 2 = 4.66 x 10“ 16 , 

k (t) = { exp( ^) )! if sin ^) > 0 
\ 0, otherwise 

for i = 3,4, with a 3 = 22.62, a A = 7.601 and uj = 43 2 00 . The constant 43200 
is 12 h measured in seconds. The initial conditions are ui(0) = 10 6 and 
2 / 2 ( 0 ) = 10 12 . 

This problem has important features like: 

• The Jacobian matrix is not a constant. 
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Figure 1: First component of the Chapman atmosphere problem 


• The diurnal effect is present. 

• The oscillations are fast. 

• The time interval used is fairly long, 0 < t < 8.64 10 5 , or 10 days. 

Let h > 0. Given different coefficients c*, 1 < i < s there is a (unique 
for h sufficiently small) polynomial of collocation q(t) of degree less than or 
equal to s such that 

q(t 0 )=yo, q(t 0 + Cih) = f(t 0 + Cih, q(t 0 + Cih)) if 1 < i < s. (27) 

The collocation methods are defined by an approximation y(t ) ~ g(t), and 
are equivalent to implicit RK methods of s stages 

S 

h = / (to + <H h, y 0 + h Y ciij kj ), 

s i=1 (28) 

2/i = y 0 + h'£b i k il 
1=1 
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Figure 2: Second component of the Chapman atmosphere problem 


for the coefficients 


a i,j — 


k = 


ret 

n 


u- Cl 


' 0 Wi C i~ Cl 


du , 


n 


u- Cl 


(29) 


du. 


' o Ci ~ Cl 


The coefficients c* play the role of the nodes of the quadrature formula, and 


the associated coefficients 6 * are analogous to the weights. From (29) we can 


find implicit RK methods called Gauss of order 2s, Radau IA and Radau 
IIA of order 2s — 1 and Lobatto IIIA of order 2s — 2. See [15] for more 
details. 

We consider the implicit fourth order Gauss method (s = 2 as collocation 
method). We approximate the associated nonlinear systems of equations 
using our Moser-Steffensen’s method as a black box. We obtain a good 
approximation as we can see in Figures [I] and [2j Note that 1/2 = [O 3 ] looks 
like a staircase with a rise at midday every day and yi = \0\ looks like a 
spike with its amplitude increases each day. 
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